Density fingerprint of giant vortices in Fermi gases near a Feshbach resonance 
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The structure of multiply quantized or giant vortex states in atomic Fermi gases across a Feshbach 
resonance is studied within the context of self-consistent Bogoliubov-de Gennes theory. The particle 
density profile of vortices with k > 1 flux quanta is calculated. Owing to k discrete branches of 
vortex core bound states, inside the core the density oscillates as a function of the distance from 
the vortex line and displays a non-monotoic dependence on the interaction strengths, in marked 
contrast to the singly quantized case, in which the density depletes monotonically. This feature, 
never reported so far, can make a direct visualization of the giant vortices in atomic Fermi gases. 
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One of the hallmarks of superfluidity of quantum fluids, 
be it fermionic or bosonic, is the appearance of quantized 
vortices. In condensed matter physics, the vortex states 
have been studied widely in various systems, ranging 
from the conventional Bardeen-Cooper-Schrieffer (BCS) 
superconductors to the rotating helium superfluids. The 
recent manipulated ultracold atomic 6 Li and 40 K gases 
emerge as a new promising test-bed for the vortex physics 
[H. Their interactions can be arbitrarily and precisely 
enhanced using a Feshbach resonance. Upon sweeping 
a magnetic field downward through the resonance, these 
Fermi systems undergo a smooth crossover from the BCS 
superfluidity to the Bose-Einstein condensation (BEC) 
of tightly bound pairs [l|. As the underlying statistics 
of systems changes from fermionic to bosonic across the 
resonance, it is interesting to ask how the properties of 
vortices evolve around the crossover. 

A singly quantized vortex in the BCS-BEC crossover 

has been discussed to a certain extent i a a a a a 

a El- The presence of strong interactions is shown to 
lead to a significant depletion of the particle density in 
the region of the vortex core 0, [H]> which was indeed 
confirmed experimentally by Zwierlein et al. for a 6 Li 
gas The properties of giant vortex states with mul- 
tiple flux quanta EH El! on the other hand, is less 
known. Generically, in a bulk system giant vortices are 
not energetically favorable and are not expected to per- 
sist if created. In the confined geometry, however, the 
situation may be different. A number of methods have 
been proposed to overcome this vortex dissociation insta- 
bility, including the use of an external repulsive pinning 
potential [l3| or a trapping potential steeper than the 
harmonic traps [14]. As a counterpart, the giant vortex 
structures have been recently produced in rapid rotat- 
ing trapped BECs They have also been observed 
in the nanoscale superconductors where the sample size 
becomes comparable to the superconducting coherence 
length £ [Hj]. 

Given all the recent advances in experimental tech- 



niques, in this paper we discuss the evolution of the gi- 
ant vortex structure from weak-to-strong coupling super- 
fluidity in trapped Fermi gases across a broad Feshbach 
resonance. In marked contrast to the singly quantized 
vortex, we find a non-trivial oscillation behavior in the 
particle density profile of giant vortices inside the core in 
the strongly interacting BCS-BEC crossover regime. The 
oscillation pattern, unique to the number of flux quanta 
at particular couplings, relates directly to the microscopic 
electronic spectrum of the local density of states (LDOS), 
which acquires an intriguing structure owing to the mul- 
tiple branches of vortex core bound states, the so-called 
Caroli-de Gennes-Matricon (CdGM) states In this 
respect, it provides a density fingerprint for giant vortices 
in the neutral Fermi gases. Towards the deep BEC limit, 
these oscillation patterns cease to exist, and finally the 
density profile returns back to that of a BEC. Our results 
are obtained by numerically solving the Bogoliubov-de 
Gennes (BdG) equations in a fully self-consistent fash- 
ion. As the strongly interacting Fermi systems can be 
found also in various fields of physics, such as the high- 
temperature superconductors and neutron stars, our re- 
sults can have implications beyond the cold atom physics. 



To be concrete, we consider a two-dimensional (2D) 
Fermi gas that can be prepared readily in a single 'pan- 
cake' trap or at the nodes of ID optical lattice poten- 
tials. It is sufficient to model the broad Feshbach reso- 
nance using a single channel Hamiltonian [l3|. We there- 
fore assume a 2D contact interaction parameterized by 
a coupling constant g. The two-body interaction prob- 
lem under this circumstance involves two length scales: 
the characteristic length in the tightly confined direc- 
tion ao and the 3D s-wave scattering length a sc . A 
peculiar 2D bound state of two atoms appears for an 
arbitrarily weak attraction [l0], with the binding en- 
ergy E a /(hu>o) = 0.915/7rexp(— \/27rao/a sc ) -C 1, where 
ujo = fi/(ma§). The bare coupling constant can then be 
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regularized via the s-wave scattering phase shift [2( 
1 
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i.e., 



(1) 



where the relative collision energy E is of the order of 
the Fermi energy Ep and drops automatically out of 
the final results. For a uniform gas at zero tempera- 
ture, the mean-field theory of the BCS-BEC crossover in 
2D admits simple analytic expressions for the order pa- 
rameter and chemical potential: A = (2EpE a ) 1 ^ 2 and 
fi = Ep — E a /2, respectively [13]. Hence, E a <C Ef 
corresponds to the weak coupling BCS limit, while in the 
opposite BEC limit of very strong attractions, E a 3> Ef. 
The crossover occurs approximately at E a ~ O.bEp [21(. 

In BdG approach the quasiparticle wave functions u v 
and v v are determined by the coupled equations [22l |. 
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where E v is the excitation energy, and the single particle 
Hamiltonian in traps is Tia = — h 2 V 2 /2m + muj 2 r 2 /2 — 
fx. As the BdG equations are invariant under the re- 
placement iijj (r) — > v* (r), v n (r) — > — w*(r), E 1 ^ — > 
— -Ejj, we restrict our calculations to E n > only. 
The order parameter A(r) and the chemical potential 
/i are determined respectively by the self-consistency 
equation A(r) = g u v (r) v* (r) [1 - 2/(J5,)] and 

the particle density n (r) = 2j^ (r)| 2 f(E n ) + 

\v v (r)| 2 [1 -/(#,,)]} so that / dm (r) = AT. Here / (x) = 
1/ {e x l ksT + l) is the Fermi distribution function and JV 
is the number of total atoms. Numerically one has to 
truncate the summation over the energy levels. In prac- 
tice, we develop a hybrid procedure by introducing a high 
energy cut-off E c , above which we use a local density 
approximation (LDA) for the high- lying modes. Thus 
the regularization prescription (U) leads naturally to an 
effective coupling constant in the self-consistency equa- 
tion A(r) = g eff (r) J2 V to v* (r) [1 - 2f(E n )], where 
53 is now restricted to E n < E c . Further expression of 
g e ff (r) and the detailed LDA contributions to the par- 
ticle density will be reported elsewhere. Below E c , we 
solve the BdG equations by taking A(r) = A(r)e~ 4K(/ ', 
where k denotes the number of vortex flux quanta. Ac- 
cordingly, we write, for the normalized modes, u v (r) = 
u nm (r) e^Vv^F and v v ( r ) = V nm (r) e^ m +^ / V2^. 
The BdG equations then decouple into different m sec- 
tors and reduce to a matrix diagonalization problem if 
one expands u nrn (r) and v nm (r) in a basis set of 2D 
harmonic oscillators. 

We have performed self-consistent computations for a 
gas with N = 1000 for k up to 10, and have set a^o — 
{h/muj) 1 / 2 and huj as the units of length and energy, re- 
spectively. In the absence of vortices, an interesting as- 
pect of the 2D mean-field theory is that the density pro- 
file is essentially independent on the interactions, though 
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Figure 1: (Color online). Particle density profiles, normalized 
by n T F = vN/ira^, for several values of the interaction 
strengths: E a = O.IEf (solid lines), E a = Q.2Ef (dashed 
lines), E a = 0.5Ef (dotted lines), and E a = 2.QEf (dash- 
dotted lines). 



S °- 5 
s 



. «0 ' O A. 










i (a) 


(b) 





"0 2 4 6 
K 



1.0 



0.5 



-4 -3 -2-10 1 2 
ln(2£ IE) 



Figure 2: (Color online). Centre particle density as a function 
of the number of flux quanta (a) and the interaction strengths 
(b). The symbols in the left panel denote different value of 
the interaction strengths: E a = O.IEf (BCS side, squares), 
E a — 0.5E F (crossover point, circles), and E a = 1.5E F (BEC 
side, triangles). The lines in the right panel show the results 
with different number of flux quanta: K = 1 (solid line), k = 2 
(dashed line), and k = 3 (dash-dotted line). 



the chemical potential is appreciably reduced. Within 
LDA we find that n (r) K=0 = (Vn/f)(1 - r 2 /r^ F )a^ 
with ttf — V / 2A 1 / 4 a/i , and fi K= o — Ep — E a /2, where 
Ep = \/~Ntvjo = ksTp. The resulting maximum value of 

1 /2 

the order parameter is Ao = (2E a /Ep) 1 Ep. We have 
chosen E c ~ 4:Ep, which is already sufficient large to 
ensure the cut-off independence of our results. The char- 
acteristic length scale of the core size of giant vortices 
is £ K ~ k£, where the coherence length ^ ~ Hvf/ttAq ~ 
kp 1 = x /l/2N- 1/4 a ho in the BCS-BEC crossover regime. 

Our main results are summarized in Figs. 1 and 2 
where we report the vortex particle density profiles and 
centre particle densities for several values of the number 
of flux quanta and the interaction strengths at nearly zero 
temperature T — 0.011>. The most unexpected feature 
of these profiles is the prominent oscillation behavior in 
the region of the vortex axis for giant vortices with k > 2, 
in sharp contrast to the monotonic depletion of the parti- 
cle density in case of a singly quantized vortex as shown in 
Fig. la. In addition, the centre particle density oscillates 
with k and displays a non-monotonic dependence on the 
interactions. The oscillations in the density profile are 
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most pronounced on the BCS side and at large number 
of flux quanta. However, they get suppressed appreciably 
with increasing the interaction strengths. Nevertheless, 
they are clearly visible around the crossover regime, and 
should be easily detected by the absorption imaging in 
experiments. For a given interaction strength, there is a 
critical value of the number of flux quanta required to ex- 
hibit the oscillations, which increases as the interaction 
increases. In the nearly BEC regime at E a = 1.5Ef, 
the oscillation occurs for n > 6 only, as displayed by 
the triangles in Fig. 2a. We thus expect that in the ex- 
treme BEC limit, these oscillation patterns in the density 
profiles of all giant vortices should vanish identically, in 
accordance with the general picture of a fully condensed 
BEC. 

The appearance of the intriguing oscillations in the 
particle density profile for giant vortices is in close con- 
nection to the multiple branches of CdGM bound states 
inside the vortex core [l3|. In the weak coupling limit, a 
simple semiclassical treatment of the CdGM states leads 
to a linear spectrum [H], 



(n + 



)E K o 



)kE k1 , 



(3) 



where to is the angular momentum, and the branch in- 
dex n may take k integrate values, i.e., from — k/2 to 
(k — l)/2, according to the index theorem established 
by Volovik for the number of anomalous branches of 
low-energy quasiparticles inside the core [23j. E K o — 
Trhv F /(2^ K ) ~ 2A /k and E Kl = h 2 /(2m£V) - 
(AI/Ef)/k 2 are the bound state level spacings [121 ]. 

To illustrate the relation between CdGM states and 
our results on the particle density profiles, we calculate 
the LDOS N(r,E), given by 2 E„[K(r)| 2 5(E - E v ) + 

li^fr)! 2 5(E + E v )}, which, when integrated over negative 
energy, gives rise to the density profiles n (r). The CdGM 
states would exhibit themselves as peaks in the LDOS. 
As the radial functions behave as u nm (r) ~ r \m+K/2\ anc j 
Vnm (r) ~ r \m-K/2\ c i ose to the origin, the quasiparticle 
probability amplitudes |w(r)| 2 and |w(r)| 2 have maxima 
at r ~ [mj /kp because of the angular momentum of the 
states [10(1 . Therefore, the principal contribution to the 
LDOS at given (r, E) arises from the bound states with 
(|m|,e nm ) = {k F r,E) 0. 

Let us first focus on the centre particle density with 
to ~ 0. Fig. 3 shows the LDOS at the vortex axis for 
different values of the interactions. On the BCS side, 
where E k q ^> E K \ (see, i.e., Fig. 3a), there are peaks lo- 
cated both below and above the Fermi surface of E = 
for K > 2, and their weights may change periodically as 
a function of k. As a result, the centre density oscil- 
lates with the number of flux quanta. With increasing 
the interactions, however, the level spacing E K i becomes 
progressively larger due to the enhancement of Ao, and 
therefore E k q < E K \ across the crossover point. Hence, 
the interaction turns to expel the bound states towards 
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Figure 3: (Color online). Local fermionic density of states at 
the vortex axis for (a) E a = O.lEp, (b) E a — 0.5E F , (c) and 
E a = 1.5E F . 
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Figure 4: (Color online). Spatial variations of the local den- 
sity of states for k — 4 and K = 5 at several values of the 
interaction strengths as labeled. 



the positive energy side. This results a sudden drop of 
the centre density at a critical coupling strength once the 
lowest bound state shifts up across E = 0, as shown in 
Fig. 2b. The value of the critical coupling increases with 
the number of flux quanta. 

We now consider the oscillations in the particle density 
profiles of giant vortices, which may be understood from 
the spatial dependence of the LDOS, as displayed in Fig. 
4 for k = 4 and k = 5. In the weak coupling limit, the k 
branch spectra of CdGM states are quasi-continuous. It 
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Figure 5: (Color online). Order parameter profiles and cur- 
rent distributions at the vortex core for several values of the 
interaction strengths: E a — O.IEf (solid lines), E a — 0.2Ef 
(dashed lines), E a = 0.5-Ef (dotted lines), and E a = 1.5Ef 
(dash-dotted lines). The currents are in units of titfVf, 
where vf = kkp/m is the Fermi velocity. 



free motion of atoms in a box of length L in z axis, we 
have also studied the 3D situation for a strongly interact- 
ing gas of N — 10 4 atoms in a cylinder with L ~ /r TF , 
and have observed qualitatively the same features. 

In conclusion, by self-consistently solving the mean- 
field Bogoliubov-de Gennes equations we have analyzed 
the structure of giant vortices in a superfluid atomic 
Fermi gas in the strongly interacting BCS-BEC crossover 
regime. The multiple branches of the CdGM bound 
states are shown to have a significant impact on the lo- 
cal density of states, and consequently lead to nontrivial 
oscillations in the giant vortex density profiles. These 
distinct oscillations, which could be visualized after ex- 
panding the cloud, can make a useful diagnosis of giant 
vortices in atomic Fermi gases. 
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and by the National Science Foundation of China under 
Grant No. NSFC-10574080 and the National Fundamen- 
tal Research Program under Grant No. 2006CB921404. 



is easy to see from Eq. ((3j) that a wedg e-shaped pattern 
of maxima in the LDOS is formed [lot. Tllj| . as shown in 
Fig. 4a. There are k rows of peaks as one moves away 
from the vortex axis, with decreasing number of peaks 
one by one due to the increase of the angular momentum 
to. Therefore, the integration over the negative energy of 
the LDOS naturally yields the oscillation behavior of the 
density profiles. However, as noted above, the increase 
of the interaction will make the CdGM states more dis- 
crete, with a larger level spacing. This destroys gradually 
the regular pattern of maxima in the LDOS and the re- 
sulting oscillations in the density profile. For a sufficient 
attraction, see, i.e., Fig. 4c, the LDOS exhausts at the 
negative energy, and therefore the density profile of giant 
vortices depletes completely inside the core, resembling 
that of an ideal BEC as expected. 

Finally, in Fig. 5 we report the order parameter 
profiles and the current circulating around the vortex 
core. Formally the current density is given by j (r) = 
(2ih/mr)^2 r] v v d i pV*f(—E r] )ip. The order parameter in- 
side the core expands as the flux quanta increases, in 
accordance with the asymptotic Ginzburg-Landau form 
A(r) ~ r K (r < £ K ), and enhances with increasing the 
strength interactions. On the other hand, the current 
density exhibits a similar oscillation behavior as the par- 
ticle density for weak interactions. These oscillations 
are attributed to the interplay between the paramag- 
netic bound states and the diamagnetic scattering states, 
which give the opposite contributions to the current, as 
discussed in Ref. [24| for a two-quantum vortex. 

We so far confine to the 2D geometry. By allowing a 
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